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Pair densities in density functional theory 


Huajie Chen* and Gero FriescckP 


Abstract 

The exact interaction energy of a many-electron system is determined by the electron 
pair density, which is not we 11-approximated in standard Kohn-Sham density functional 
models. Here we study the (complicated but well-defined) exact universal map from 
density to pair density. We survey how many common functionals, including the most 
basic version of the LDA (Dirac exchange with no correlation contribution), arise from 
particular approximations of this map. We develop an algorithm to compute the map 
numerically, and apply it to one-parameter families {ap{ax )} a >o of one-dimensional 
homogeneous and inhomogeneous single-particle densities. We observe that the pair 
density develops remarkable multiscale patterns which strongly depend on both the 
particle number and the “width” aT 1 of the single-particle density. The simulation 
results are confirmed by rigorous asymptotic results in the limiting regimes a » 1 and 
a << 1. For one-dimensional homogeneous systems, we show that the whole spectrum 
of patterns is reproduced surprisingly well by a simple asymptotics-based ansatz which 
slowly smoothens out the ‘strictly correlated’ a = 0 pair density while slowly turning 
on the a = oo ‘exchange’ terms as a increases. Our findings lend theoretical support to 
the celebrated semi-empirical idea to mix in a fractional amount of exchange, albeit 
not to assuming the mixing to be additive and taking the fraction to be a system- 
independent constant. 


1 Introduction 

Density functional theory (DFT) |23| [241, 3J4 provides the most widely used models for com¬ 
puting ground state electronic energies and densities in chemistry, materials science, biology, 
and nanosciences. The success of DFT lies in the use of exchange-correlation functionals 
that model the intricate many-body interaction energy by explicit expressions in terms of 
the one-body density or the one-body Kohn-Sham orbitals. Although currently available ap¬ 
proximations, such as B3LYP sms] or PBE [34], perform remarkably well for a wide range 
of systems, the DFT models exhibit well-known failures when strong correlation effects are 
present, as arising for example in the breaking of chemical bonds [9]. Therefore, finding 
an accurate single-particle formalism that remains reliable in strongly correlated regimes 
remains a major challenge. 

In this paper we shed new light on this challenge by studying the exact map p D p 2 from 
single-particle density to pair density whose existence is assured by abstract DFT. The exact 
interaction energy is obtained by integrating the pair density against the Coulomb repulsion 
potential (see <E1 below), so any approximation p i— p 2 yields an approximate interac¬ 
tion energy functional. We take the view, first advocated by Gunnarsson and Lundqvist 
pi] , that the exact density-to-pair-density map p i-T p 2 is a better starting point to under¬ 
stand or design model interaction energy functionals than the commonly used density-to- 
interaction-energy map p V ee [p ]. This is because the pair density, a function on two-body 
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configuration space, encodes a wealth of physically and mathematially interesting informa¬ 
tion about a many-body quantum system which is “averaged out” in the interaction energy, 
a mere number. In particular, comparing exact and approximate pair densities does not just 
yield a total interaction energy error, but also reveals where in two-body configuration space 
the error is localized. 


The main results in our paper are careful simulations of the exact density-to-pair-density 
map for typical one-dimensional model systems, with different electron numbers and different 
density profiles varying from “concentrated” to “dilute”. The algorithm we develop for this 
purpose allows one to deal with the infinite-dimensional, nonlinear constraint of fixed single¬ 
particle density which appears in the definition of the map. To obtain a simple form of this 
constraint after discretization, we use a finite-element basis from computational mathemat¬ 
ics instead of the usual basis sets of quantum chemistry. We observe remarkable multi-scale 
patterns in the pair density which strongly depend on both the particle number and the 
“width” of the density profile. See e.g. Figure |(h9| in Section [6^2) These patterns are not 
accurately captured (except in extreme regimes) by any of the currently used DFT models, 
but in our view constitute fundamental low-dimensional manifestations of exact DFT. We 
thus hope that our simulations, despite the limitation to one-dimensional model densities, 
offer exciting glimpses of possible future DFT models. 


In the remainder of this Introduction we informally discuss the definition of the exact density- 
to-pair density map, display two instructive extreme and opposite approximations in the 
DFT literature, and explain how our simulations seamlessly connect all three. The exact 
density-to-pair-density map is constructed as follows: 


p i—>• T p2 (1.1) 

where p T is the map obtained by Levy-Lieb constrained search [26] [29], i.e. T is 
the Welectron wavefunction which minimizes kinetic plus potential energy, T[T] + V ee [T], 
subject to the constraint that T has single-particle density p (see Section 2 for notation, 
function spaces, and further explanation), and the pair density associated to any Welectron 
wavefunction T is 

|T(x,cri,7/,cr 2 ,X3,cr 3 , • • • , x N , (J N )\ 2 dx 3 • • • dx N . (1.2) 



Here (aq, cq) G M 3 xZ 2 are space-spin coordinates for the i th electron. Minimizing T’s always 
exist [29 , and the complication of possible non-uniqueness is discussed in Section [3j The 
electron-electron interaction energy is a simple explicit functional of the pair density, 

V ee m= f f^- dxdy. (1-3) 

J R6 \x-y\ 


Hence any approximate expression of the pair density in terms of the single-particle density 
gives, by substitution into (1.3), an approximate interaction energy functional. Numerous 
functionals have been formulated in this way eu la mum ins]. Standard DFT models start 
from a statistical independence ansatz 


P2(x,y) = -p(x)p(y), 


(1.4) 


and include all the many-body effects in a correcting exchange-correlation energy functional. 
Opposite to this uncorrelated ansatz, there is the more recent strictly correlated electrons 
(SCE) model [37, 38, [40] . which is attracting attention in the mathematics literature [TT 
mmm due to its connection with optimal transportation theory and which arises from 
neglecting the kinetic energy in the constrained search in The corresponding ansatz 

for the pair density is 


P2(x,y) = T7E / P(z)5(x-Ti(z))6(y-Tj(z))dz 
2N j M3 


(1.5) 
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with Ti : M 3 -A M 3 , i = 1, • • • } N being certain optimal transport maps, see Section [ 4 ] for 
more details. 

The formalisms <0 and ( |1.5| ) give “extreme” pair densities. True pair densities, unlike 


(1.4), are expected to localize in certain regions due to shell structure or ionicity avoidance; 


but (1.5) emphasizes this localization too much and misses its quantum features. Most of 
the practically interesting models, such as the local density approximation (see Section [4|, 
lie “inbetween” these two extreme distributions. But what kind of “interpolation” is the 
right one, and captures true pair densities 0? 

A marvellous tool to approach this question is density scaling , introduced in the context 
of exact DFT by Levy and Perdew m, which is closely related to the adiabatic connection 
utilized in many DFT studies (see e.g. 2F 27, 38, 4CK m • Alongside a given density 
p : -A M, consider - as we shall in our simulations - its re-scalings 


(D a p)(x) = a d p(ax ), a E (0, 00 ). 


( 1 . 6 ) 


The parameter a seamlessely rescales a dilute system (a « 1) into a concentrated one 
(a » 1). But the associated pair densities do not just change by a rescaling, that is to 
say p 2 [D a p\ 7 ^ D a p 2 [p\ (where D a acts on pair densities as (D a p 2 )(x,y) = a 2d p 2 (ax, ay)). 
Instead, it follows from the arguments in [27 that the following diagram commutes: 


scale 


min aT+Vee 


D a p 

min T+Ve e 


(1.7) 


7—v r t ~\ "I SCctl.6 back r ta 1 

D a -ip 2 [D a p\ ^ - P2 [D a p\ 

Here ‘min’ means find the minimizing wavefunction under the constraint of the given one- 
body density and take the resulting pair density. See Proposition |5.1| below. Thus the scaling 
parameter a in (1.6) acts as a coupling constant in the one-parameter family of variational 
problems on the left which govern the pair density. This family “adiabatically”, i.e. while 
keeping the density fixed, connects the problem of minimizing just V ee (a = 0) via T + V ee 
(a = 1) to minimizing just T (a = 00 ) FI Nontrivial but well known formal asymptotics for 
the minimizing wavefunction for a -A 0 (37) and a -A oc (see e.g. HD together with formula 
0 then suggests the following: the true pair density is asymptotic to that of the SCE 
state, eq. (1.5), as a -A 0 S3, and to that of the Slater determinant of the Kohn-Sham 


orbitals as a -A 00 . The latter reduces to 0 when N = 2 or when the particles are bosons, 
and in addition contains ‘exact exchange’ (see Section |4| for higher N. See Sections HI IT- 
more details and rigorous proofs in special cases. 

Now back to our central question: which “interpolation” between the extreme pair densi¬ 
ties (|1.4|) and (|1.5|) is right? Our numerical results for typical families (1.6) of one-dimensional 


densities with different particle numbers show that many different interpolations are right. 
See, e.g., Figure [679] The true pair densities form a two-parameter family which strongly de¬ 
pend on both the particle number N and the scaling parameter a. At fixed N they steadily 
“cross over” from ( |1.4|) (plus exact exchange when AT > 2) to (1.5). The impractical, highly 
implicit definition ( |1.2| ) is able to pick out the right parameter values from the density, but 
simple explicit formulae will not. In the very special case of homogeneous systems in one 


dimension (see Figure 6.8) we design an ansatz which does. The idea is to simultaneously 
smoothen out the pair densities from the strongly interacting limit and fading out the ex¬ 
change terms from the weakly interacting limit. But the correct smoothing lengthscale and 
the correct fraction of exchange keep changing with N and a. See Table [3j Our simulations 
to some extent support the celebrated idea pFl underlying the functional B3LYP to mix in a 


Andreas Savin suggested to us the name two-sided adiabatic connection because it combines the classical 
connection to T (in our parametrization, 1/a E [0,1]) with the more recent one to V ee |40| («e [ 0 , 1 ]). 

2 According to a recent article in Nature (29.10.2014), one of the Top Ten most highly cited scientific 
papers of all time, and the most highly cited one written after 1990. 
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fraction of exact exchange. But they show that the right fraction, taken to be 0.2 in B3LYP 
[5], is in fact not constant. At present we have no proposal how the right fraction could be 
adaptively picked out in realistic (inhomogeneous, 3D) simulations. 

The remainder of this paper is arranged as follows. In the next two sections, we recall 
basic aspects of DFT and give the precise definition of the universal density-to-pair-density 
map. We then show in Section [4] how some common DFT functionals arise from approxima¬ 
tions of this map. Section [6] describes our numerical simulations of the true pair densities 
of homogeneous and inhomogeneous one dimensional systems for both bosons and fermions. 
In Section [7] we present rigorous asymptotic results which confirm the numerical findings. 
In Section [8| we propose an ansatz for approximating pair densities of one dimensional 
homogeneous electron systems. Finally, in Section [9] we give conclusions and some future 
perspectives. 

2 Density functional theory 

Here we recall the basic functionals of DFT which will be needed in the following. A standard 
reference is [33]. Readers familiar with DFT might want to skip this section. We consider 
a general system of N nonrelativistic electrons in R d under the influence of an external 
potential v ext : R d M and a repulsive pair potential v ee : R d M. Prototypically, for 
real physical systems, 


d = 3, v ee (x - y) = - —-— T (x, y eR d ), 

\x - y | 

and v ext is the electrostatic potential generated by M nuclei of charges Zi,.., Zm > 0 located 
at positions i?i,.., Rm E M 3 , 


M 


(* eRS >- 

The quantum mechanical ground state energy of the system is given by 


E 0 = ^inf (r[tf] + V ee [*} + V[*]), 


veA N 

where An is the following class of admissible wavefunctions 

An = {4/ e L 2 ((R d x Z> 2 ) N ), VT G T 2 , T antisymmetric, ||T || L 2 = l} , 
Z >2 = an d T, We, V are the following functionals: 

i r N 

r [^] = o / T V \V x N(xi,°i, ••• ,X N , (Jiv)| 2 ^i ■■■dx N 

2J (® d ) N * 1 ,-,*Nez 2 i =i 


(kinetic energy), 


V £ 


( 2 . 1 ) 


( 2 . 2 ) 


(2.3) 


E E v ee (xi - x j )\^(x 1 ,a 1 , ■ ■ ■ ,x N ,a N )\ 2 dxi ■ --dx N (2.4) 

,<TtvGZ2 l<i<j<N 


m = / 

Jq e 

(electron-electron interaction energy), and 

AT 

vm = [ 

l 

(external potential energy). 


E Et>e*(si)|tt(*i,qi, ' ~' ,x N ,a N )\ 2 dxi ■ --dx N (2.5) 

<r u -,* N ez 2 i= 1 
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A central result of DFT going back to Hohenberg and Kohn is the following. We state 
the result here in the form discovered by M.Levy [26] and made rigorous by |29] . The 
quantum mechanical ground state energy (2.1) can be recovered exactly by minimizing a 
certain density functional, 


E 0 = inf (Fhk[p\ + [ 
pen N \ j R 


VextP 


where 


R d 


F HK [p]= min {TW+Veel*]}. 


( 2 . 6 ) 

(2.7) 


Here T i—)> p means that T has single-particle density p, i.e. 

p(x) = N E \V(x,<j 1 ,X2,(T2,--- ,x N ,(J N )\ 2 dx 2 ---dx N , (2.8) 

J(RS)N-t _ 


cr ,a N eZ 2 


and 7Zn is the space of densities arising via (2.8) from wavefunctions T G An- Note that 
the space 1Z at of densities is known explicitly: by a result of Lieb [29] . 

H N = :R d ^R\p>0, J p = iV | . (2.9) 

We also note that Fhk is a universal functional of p, in the sense that it does not depend on 
the external potential v ext . Minimizers in (2.7) always exist provided p G 7Zn m- 

The complexity of the DFT model (2.6) lies in that no tractable expression for Fhk is 
known that could be used in numerical simulations. In practice, Fhk[p] is approximated by 
the sum of a kinetic part and an interaction part, 


Fkk\p] ~ T[p\ + V ee [p], 

leading to an approximate expression for the ground state energy, 


Eq ~ Eq — 


inf (t[p\ + V ee [p\ + f v ext p V 

PGTCn \ JR d ' 


( 2 . 10 ) 


( 2 . 11 ) 


Many clever and useful approximate functionals T and V ee have been proposed and utilized 
to simulate a wide range of systems (see e.g. [33;, 0]). Particularly fruitful has been the 
idea of Kohn and Sham [24] to construct a kinetic energy functional T with the help of 
single-particle orbitals of a non-interacting reference system: 


T[p] = T K s [p\ = min 


[lEj |V^| 2 , 


x ^ 2 ), 


N 


— ^ij •> 


E E i‘M x ’ <j )i 2 = p(x) 

i =1 ctGZ 2 > 


( 2 . 12 ) 


where for any function / = f(x, a) of (. x , a) G W l x Z2 we use the abbreviation f f = 
5D<rez 2 fu d f( x ’ a ) dx - ^ is easy to see that T is the same as the functional obtained by 
omitting V ee in (|2.7| and restricting the minimization to Slater determinants 


^(xi,a u ..,x N ,a N ) = 


Vm 




<Pn(x 1 , 01 ) 


<Pi(xn,ct N ) ■■■ <pn(xn,<?n) 


(2.13) 


Minimizers in (2.12) always exist provided p G 7Zn m- 
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Using (2.12) as the kinetic energy in (2.10), as is done in almost all simulations to date, and 
the orbitals = (</>i, • • • , n) as the basic variable, the ground state energy of the system 
becomes 


E 0 nE 0 = mf\l V [ |V^| 2 + [ 

{ 2 i=l J • / K° 


with the single-particle density 


VextP& T h^ e [p<|)], (f)i G H 


'(m 3 ), /; 


~ $ij 


(2.14) 


N 

p ^{ x ) =J2J2 i(/>i(x, s )i 2 . 

^—1 

The remaining problem, and the one of interest to us, is to design accurate approximations 
for y ee [/>$]. 


3 Universal density-to-pair-density map 


Our starting point for looking at interaction energy functionals will be the universal, exact 
density to pair density map delivered by abstract DFT. Following Levy [26] this map is 
defined as follows. Recall from Ql.2| ) that p% denotes the pair density of the wavefunction T. 


Definition (Universal density to pair density map) For any one-body density p of an N- 
electron system, that is to say for any p belonging to the class 1Zjy in (2.9), 


P 2 [p\ = {p 2 | T G *4at is a minimizer of T + V ee subject to T p}. (3.1) 


Just like the map p i-A Fhk [p\ , the map p \-> p 2 [p\ is universal, i.e. independent of the 
external potential. The above definition requires, and it was proved mathematially by Lieb 
[29] . that a minimizing T exists. Note however that the minimizer may not be unique. 
Hence the map is possibly multi-valued, that is to say p 2 [p\ is possibly a set of pair densities 
rather than a single pair density. Simple explicit examples of nonuniqueness in the case when 
T + V ee is replaced by T are given in Section [7| 


The physical significance of p 2 [p\ comes from the following direct consequence of formulae 
(2.1), (2.6): if T is any exact quantum mechanical ground state, i.e. a minimizer of the right 
hand side of (2.1) for some external potential v ext , and T has one-body density p, then p 2 [p\ 
is the exact pair density of T, and the functional 


Vee[p}'- v ee {x-y)p 2 [p\(x,y)dxdy 

J R d 


agrees with the exact interaction energy Ueel^] from (2.4). 


(3-2) 


4 Approximate density-to-pair-density maps 


It is obvious that substituting any approximation p 2 [p\ of the density to pair density map 
P 2 [p] into (3.2) yields an approximate interaction energy functional V ee [p\. Conversely, we 
now show that many basic approximate functionals used in practice can be derived in this 
way. In some cases, such as the bare Hartree functional or ‘exact exchange’ (Examples 1 and 
3), this is trivial. For the LDA in its most basic form (Dirac exchange with no correlation 
contribution, Example 2) it is not, and we are not aware that an exact equivalence to a pair 
density model for any inhomogeneous density as given below has been stated previously, even 
though good approximate pair density formulations are well known ED- For interesting work 
relating advanced DFT functionals to pair density approximations we refer to |2j ESI [35] . 
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Example 1. (statistical independence) The simplest idea is to assume statistical indepen¬ 
dence, 


P 2 [p]{x,y) = -p(x)p(y). 


(4.3) 


Substituting this density to pair density map into the right-hand side of (3.2) leads to the 
Hartree functional 

Vee[p] = \ [ p( f} p ^ dxdy. (4.4) 

2 Jr e \x - y | 


While never used on its own, together with some correcting exchange-correlation functional 
E xc [p\ it is contained in virtually all DFT models, including state of the art ones like B3LYP 
|25l 3 or PBE [34] . 


Example 2. (Local density approximation with Dirac exchange) For the free (i.e., nonin¬ 
teracting) electron gas, the pair density can be determined explicitly (see e.g. [33] and, for a 
mathematical account, [15j). In this case the single-particle density is a constant, p(x) = p, 
and the pair density is 


P 2 (x, V ) = ^P 2 ~ ^p 2 h 2 ((3n 2 p) 1/3 \x - y\ ), (4.5) 

where h(s) = 3(sins — scoss)/^ 3 . We claim that the inhomogeneous version 

P 2 [p](x,y) = ^p(x)p(y) - ^p(x) 2 h 2 ((3Tr 2 p(x)) 1/3 \x - y\) - ^p{y) 2 h 2 {(3n 2 p(y)) 1/3 \x - y\) 

(4.6) 

yields the interaction energy 

Vee[p]~\ [ dxdy - c x f p(x) 4/3 dx (4.7) 

2 J R6 \x - y I J R 3 

with constant c x = l(f) 1 ^ 3 - This can be seen as follows. For each of the non-mean-field 
terms, just integrate out the variable not contained in the argument of ft, e.g., using spherical 
polar coordinates for y centered at x and abbreviating kp(x) = (3 tt 2 p(x)) 1 / 3 , 

/ h 2 (k F (x)\x-y\)dy = 4ir h 2 (k F (x)r)r 2 dr = - - - — T - / h 2 {r')r' 2 dr ', 

J R3 Jo (37 T Z p(x)) Z / 6 J 0 


and determine the remaining one-dimensional integal as in the discussion of the homogeneous 


case in [33j[l5|. Eq. (4.7) is the simplest of the local density approximations (LDA) j ~24l [33l 
eh. The second term of ( |4.7| ) is the celebrated Dirac exchange functional na. We remark 
that from Dirac’s original derivation it is not clear how to relate this functional to the pair 
density as he used a semiclassical limit argument for the (one-body) energy density per unit 
volume. 


Strange as the model (4.6) for the pair density may look, it provides a precise way to state 


what the LDA really does: the pair density is assumed to be independent at long range (note 
that ft(r) goes to zero as r gets large), while at short range it contains an “exchange hole” 
0of fixed shape coming from free electron gas theory whose diameter is of order p(x) -1 / 3 . 

Hartree-Fock-like, exchange [3], one 


i.e. 


Example 3. (exact exchange) To obtain “exact” 
takes 

h\p\(v,y)~pt(x,v), (4-8) 

where T is the Slater determinant ( |2.13| ) composed of the (p-dependent) minimizing orbitals 
Pi,..., (fN h 1 the definition of the Kohn-Sham kinetic energy functional (2.12). A more 


3 see e.g. m for more information about this semi-empirical notion 
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explicit expression for p 2 is obtained by using the well known expression for the pair density 
of a Slater determinant (see e.g. my 


h[p](x, y ) = \p{x)p{y) - \t{x, y) with t(x, y) = Y^ 


N 




i=1 


(4.9) 


Thus, just as in Example 2 the pair density naturally decomposes into a statistically inde¬ 
pendent term plus an exchange hole, but here the shape of the hole is no longer fixed but 
adapts itself to the density at hand. Expression (|4.9|) results in the interaction energy 


Veelp] = \ / 

Z JR 6 


p(x) p(y) - r(x,y) 
\x-y\ 


dx dy. 


(4.10) 


The correction to (4.4) is known as exact exchange. Note that the resulting ground state 
energy (2.11) is not quite the Hartree-Fock energy. This is because the orbitals are only 


determined via minimization of kinetic energy, rather than self-consistently accounting also 
for exchange. However, if one treats the orbitals 4> = (0i, • • • , 0 tv) as the basic variable, views 
the right hand side of (4.9) as an orbitals-to-pair-density map p2[4>](x, p), and substitutes 
into (3.2) and ( |2.14| ) one obtains precisely the Hartree-Fock energy. 


Example 4. (Hybrid models) If we take some convex combination of (4.6) and (4.9), the 
interaction energy begins to resemble, up to certain further corrections, state of the art hybrid 
functionals such as B3LYP 13113111], which are widely used in contemporary computations. 


Example 5. (strictly correlated electrons) A more recent construction is the SCE (strictly 
correlated electrons) functional |37l 38] SO] 


Vee[p\ = = 


inf [ 


p(z) 


E 




7 dz , 


(4.11) 


with the infimum taken over maps Ti,..., T/v from M 3 to M 3 which satisfy T\(x) = x and 
which preserve p, that is to say 


/ p = / p for all measurable sets A C 

JA Jt,(A) 


/A JTi(A) 

This corresponds to the following density-to-pair-density map which we call pf 1CE [p]: 

fc\p\{x,y) = P 2 CE [p\( x ,y) = AE / p(z)S(x-T i (z))5(y-'T j (z))dz (4.12) 

2N M* 

with the Ti being minimizing maps. The physical meaning of the Ti is that the position of 
one electron (at x = Ti{x)) fixes the positions of all the other N — 1 electrons (at Ti{x) with 
2 < i < N). Mathematically, the variational problem in ( 4.1 1| ) is a multi-marginal optimal 
transport problem. Minimizers are known to exist when N = 2 nan] and d = 1 m- It is 
believed (and has been proved mathematically for N = 2 [12]) that Vf e CE agrees with the 
lowest expectation of Coulomb repulsion energy with a given single-particle density p, 


V 


SCE 


[p\ = inf V ee m. 
LPJ veA N ,*^p 


(4.13) 


To derive (4.11) from (4.13), one notes that the infimum in ( |4. 13 ) is not attained in any rea¬ 
sonable wavefunction class such as (2.2) or {T G L 2 ((M 3 xZ 2 )' /v ) : T antisymmetric, ||T || L 2 = 


1}. Therefore, one needs to augment the admissible 7V-body densities pw = J2 Sl 
in (4.13) from integrable functions to probability measures, i.e. considers 


,SjSfE^2 




mm 

PN 1 ^P J R 3 N 


j 

Jr 3 


E- 

i<j 


1 


7 dp AT, 


(4.14) 
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and makes the ansatz [37. [38] 


Pn(xi, • • • > %n) = 


M?lf n*( 


Xi - T v{ i)(z))dz, 


(4.15) 


where the sum runs over all permutations V of {1, ..,7V}. Note that (4.12) is obtained by 
integrating out all but two electron coordinates from this pN- The ansatz (4.15), which 
reduces the high-dimensional problem (4.14) to a computationally feasible one, was later un¬ 
derstood mm as an instance of the mathematical belief that “Kantorovich equals Monge”, 
i.e. that optimal Kantorovich transportation plans in are induced by Monge maps for well be¬ 
haved marginal densities p (see for pioneering results and [42] for a comprehensive 

survey). 

Examples 1 to 4 are based on a non-interacting picture and treat many-body effects as 
corrections. Despite their great successes, these models exhibit known failures for strongly 
interacting systems J9l. By comparison, Example 5 takes the strongly interacting limit and 
has been proved to be good at simulating some strongly correlated model systems (e.g. 
mi [ 32 ]), but severely underestimates the true ground state energy in standard regimes (see 
e.g. the dissociation curve of the hydrogen dimer calculated in [8]). It is therefore of great 
interest to enquire as to the structure and behaviour of the true pair densities p 2 [p \. 


Density scaling, adiabatic connection, formal asymp¬ 
totics 


In order to naturally access pair densities in different correlation regimes without changing 
the “shape” of the one-body density, we will from now on look at one-parameter families 
of one-body densities obtained by rescaling a fixed reference density p : R d —>> M (see eq. 


(1.6)). The associated pair densities do notjust change by a rescaling (see the Introduction). 
This reflects the physical phenomenon that electron correlation in dilute systems (a « 1) 
is completely different from electron correlation in high-density systems (a » 1). The 
governing variational principle for the resulting constrained-search wavefunction in was 
found by Levy and Perdew m ■ As a straightforward corollary of their analysis we obtain 
the behaviour of the density-to-pair-density map under density scaling: 

Proposition 5.1. (Density scaling) Let a > 0 and let p be any single-particle density on 


M , i.e. any function belonging to the class 7 Zn- Then the diagram (1.7) commutes. In other 
words, if p2,ct[p\ denotes the density-to-pair-density map along the adiabatic connection (left 
arrow in the diagram), that is to say 

p 2 ,a[p\ •= {p 2 | 47 is a minimizer of aT + V ee on An s/to 47 \-> (5-1) 


and p 2 [p\ is the original map (3.1), then 


D a -ip 2 [D a p\ = P2,a\p]- 


(5.2) 


Proof. For convenience of the reader we include the simple proof. For any a > 0 and any 
47 £ An , 47 p, we can rescale 47 by 

4' a (a;i, • • • ,x N ) = a dN/2 ^(ax!, ■ ■ ■ ,ax N ). 

We have that 47 a belongs to An and has one-body density D a p. Moreover 

T[4g = a 2 T[<H} and V ee [* a ] = aV ee [V}. 


4 “s/to” means “subject to” throughout this paper. 
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It follows that Tq, is a minimizer of T + V ee subject to D a p if and only if ^ is a 

minimizer of olT + V ee subject to T p. By definition, the pair densities of the minimizing 
T a ’s yield the set p 2 [D a p ], whereas the pair densities of the associated \Ks give the set 
P2,a[p\- 0 


From now on, instead of considering the scaled densities (1.6) and applying the original 


density-to-pair-density map, it is more convenient for us to fix a reference single-particle 


density and vary the coupling constant a in the constrained-search problem in (5.1) (i.e., in 
the adiabatic connection) to investigate systems in different correlation regimes. 

We note that definition (5.1) stays unchanged under multiplying aT + V ee by a positive 
constant, so one might as well use T + a~ 1 V ee . In particular, one has a well-defined density- 
to-pair-density map at a = oo: 

p2,oo[p\ = {p2 : ^ £ A/v, 4 / is a minimizer of T on An s/to 4/ p}. 

It is considered well-established in the physics literature (see e.g. [39] ) that the minimizing 
wavefunction in (5.1) has the following asymptotic behaviour: 

(5.3) 


and 


4> « Slater determinant of the KS orbitals from Ex. 3, Section [4] (o >> 1) 
TV-point density of the SCE state, eq. ( 4.15| ) (a « 1). 


£ i*i 

si ,..,sat 6^2 

Taking pair densities leads to 


P2,, 


(4.9), a » 1, 


(4.12), a « 1. 


(5.4) 


(5.5) 


Complete mathematical proofs are not available for general p. It is not clear in which sense 
to measure convergence, nor what happens if the ground state is degenerate. In fact, even 
much more basic things such as existence of optimal maps or continuity of the HK functional 
have not been proved. The rigorous analysis of ID examples in Section [7] shows that things 
are not quite as simple as one might intuitively expect. For instance, in case of orbital 


degeneracies the assertion (5.3) can be true for some choices of minimizing KS orbitals but 
not for others. 

At least for N = 2 or in the case of bosons we can offer a general result. For bosons, the 
set An of antisymmetric wavefunctions has to be replaced by 

Bn = G L 2 (R d ' N ), VT e T 2 , T symmetric, 114/|| L 2 ( R d-iV) = 1}. 


(5.6) 

Proposition 5.2. Let p be any single-particle density of an N-particle system, i.e. p G 7Zn- 
If N = 2, or if the particles are bosons, then the independent pair density 


\{}~ f)p( x )p(v) 


(5.7) 


belongs to the set p 2 ,oo[p\- 

Proof. We claim that the product wave function T(xi,.. ,xn ) = n^li \/p( x i)/N, which has 
pair density (5.7), is a minimizer of T on Bn subject to the constraint T t p. To see this, 
consider a general T E Bn with T i—)► p, and estimate 


T[T] 


1 x A f 

T \V Xi 'I'\ 2 dxi.. .dx N 

1 Jm dN 



V x N\ 2 dxi... dx n 


N I" |/ Rti (iv-i) Re(^V a , 1 ^)da;2 ..-^jv| 2 ^ 

2 Jmd f Rd( .N-i)\y\ 2 dx 2 ...dx N " ' 

n r |Vs,/ Rd(J v-i) \^?dx2 ...dxjvl 2 ^^ = i r | YPpp dx 

8 J R d p{xi) 1 8 J K d p{x) 
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On the other hand, by an elementary calculation, T[\F] is equal to the expression in the last 
line. For fermions with N = 2, analogous arguments show that the Slater determinant with 
orbitals yjp(x)/N t (s), \Jp(x)/N ( s ) is a minimizer. □ 


6 Numerical investigations of the pair densities 


We now turn to the intermediate regime where a lies somewhere inbetween zero and infinity, 


and investigate numerically how the crossover between the limit behaviour (5.4) and (5.3) 
occurs. To this end we compute, for simple reference densities p, the whole one-dimensional 
family of pair densities p2,a[p] (a €= (0, oo)) along the adiabatic connection. Recall that each 
P 2 ,a arises, up to a re-scaling, as a true pair density (see ( |5.2| )). 

Due to the nontrivial (infinite-dimensional, nonlinear) constraint 4/ p and the need to 
resolve N-e lectron wavefunctions, we limit ourselves here for simplicity to one-dimensional 
reference densities p and particle numbers N = 2,3,4. We hope that our results are never¬ 
theless of some physical and chemical interest. 

Note that the one dimensional Coulomb repulsion can not be described by l/|x| since 
the latter function is not integrable near 0. We therefore use an effective potential c( \x\) 
which is obtained by integrating the Coulomb repulsion in M 3 in a thin wire over the lateral 
degrees of freedom [5]. Explicitly, 


(^ A 
c(r) = ^" ex P 



where b is a constant and erfc is the complementary error function. We set b = 0.1 in our 
simulations (see Figure |6T| ). 

Let Q = [—L, L] (with L = 5.0 in the simulations) and let N be the particle number. We 
consider two typical systems on Vt (see Figure |R2| : a homogeneous density with periodic 
boundary condition 

p(x) = ( 6 - 1 ) 

and a smoothly varying density with zero Dirichlet boundary condition 


p(x) = 


N 

2L 



( 6 . 2 ) 


Both of these two single-particle densities belong to space (2.9) (with Ir replaced by fi). 




-5 0 5 


x 


Figure 6.1: The one dimensional effective 
Coulomb potential. 


Figure 6.2: The homogeneous and inhomoge¬ 
neous electron density (6.1) and (6.2). 


For later purposes, we calculate the optimal transport maps by using the formulae in [40] 
(which were recently justified rigorously in [12j for TV = 2 and in HD for general AT, and are 
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Figure 6.3: The optimal transport maps of the density p(x) = 



described in Theorem |7.1| below) and present them in Figure [6~3] and 6.4 for the two systems 
with 2, 3, and 4 particles. 

Moreover, in one dimension it is known rigorously m that the maps Ti,T/v are cyclic, 
that is to say 

T 2 o ... o T 2 = id, T 2 o ... o T 2 = Tj (j = 2,.., N). 


n times 


j-1 times 


For related insights see D2I- This allows to simplify formula ( |4.12| ) for the a = 0 pair 
density p 2 CE • Namely, a change of variables shows that in this case the sum over j in 
(4.12) is independent of i. This together with the fact that the normalized line element (one¬ 
dimensional Hausdorff measure) ds on the one-dimensional curve graph Tj={(x J Tj(x)) : 
x E M} is given by 


ds = ^1 + Tj(x) 2 dx, 


which yields the expression 


SCE 
P 2 


1 N 

\p]( x ,y) = o Tl 


p(x) 


y/l+T<(x)* 


ds 


y —T? ( x ) 


(6.3) 


This remarkable formula shows that the maps Ti , and hence the full 7V-body SCE density, 
can be explicitly read off from the SCE pair density! 

To obtain the true pair densities of our two typical systems for finite coupling constant 
we need to simulate the constrained-search problem in (5.1). In our case this problem is 
given, for a one-dimensional single-particle density po> by 


Minimize 


E 


<Ti,..,C r 7ve^2 



N dV 2 


dxi 


+ y \^\ 2 c(\xi - Xj\) \ dxi-■ ■ dxN 
1 <i<j<N J 

s/to 'J e An, vP >-->■ po • (6.4) 
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Here for the inhomogeneous density (6.2) An is the standard wavefunction class (2.2) with 
replaced by Q = [—L, L], 


An = VT G L 2 (([L,L\ x Z 2 ) iV ), 4/ antisymmetric, ||4/|| L 2 = l|. 


(6.5) 


Since p is zero at d=L, the constraint T i->> p automatically implies Dirichlet zero boundary 
conditions T =0. For the homogeneous density (6.1), we use periodic wavefunctions 


Xi =±L 


A n = {'J : V'L e L 2 (([-L,L] x Z 2 ) n ), 


T antisymmetric, T 


= T 


Xi=L 


x; = — L 


for all i, 11 ^||l 2 = l}. 


( 6 . 6 ) 


Recall that the constraint T i—>• po means that integrating |T| 2 over all but one electron 
positions and summing over all spins gives the single-particle density po. By the symmetry 
of |T| 2 , one can leave any of the electron coordinate not to be integrated. Therefore, the 
associated Lagrange function of (6.4) is, abbreviating Zi = {pc^ <Ji) G QxZ 2 and dxi = 

fnx z 2 ^ Zi > 


jC(9,X lt --' ,X n)= [ 
J(i 


(nxz 2 ) N 


f Z \^\ 2 c(\Xi- Xj \)] d Zl ---dz N 

z i=1 1 l<i<j<N 


+ Xi(x) ( po(x) N / \'&(x,a,Z2, ■ ■ ■ ,z N )\ 2 dz 2 ■ ■ ■ dz N ) dx 

Jn \ aeZ2 J(nxz 2 )"-i ) 


+ • • • + [ ^at(x) ( p 0 (x) — n f 

Jn \ aeZ2 J(ftxz 2 y 


|\E , (;zi, • • • , zn— i,x, a)\ 2 dzi ■ ■ ■ dzN -1 ) dx 


with the Lagrange multipliers Ai (x), • • • , A.y (x). By the symmetry of '1' | 2 . we have Ai(x) = 
••• = A,y(;e) := A( x ). Therefore minimizers of ( |6.1[ ) satisfy the following Euler-Lagrange 
equation 


-^A+ y] c(\xi - xj\) I • • • ,x N ,a N ) 


1 <i<j<N 


N 


= y^A (Xi) ^(xi,cri, • • • ,x N ,a N ) 


(6.7) 


, ’L e A n , ^ po 


\i=l 


Formally, the Lagrange multiplier X(x) equals the functional derivative of the Hohenberg- 
Kohn functional (6.4) with respect to electron density, and —X(x) equals the external po¬ 
tential vq for which po is the ground state of the system, 


A = 


dF a [p\ 


dp 


= -vq 


P=P o 


with F a [p\ := mim^^^^p {<xT[T] + V ee [&]}. Therefore, by using the Euler-Lagrange 
equation J6.7| ) we implicitly require that the density po can be generated by some external 
potential jjThis is called ‘Vrepresentability” ED, and the conditions for such densities 
are not known in general. In particular, we do not know rigorously whether the single¬ 


particle densities given by (6.1) and (6.2) are f-representable. Nevertheless, after numerical 


5 We thank Eric Cances for this remark. 
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discretization it can easily be shown that the Lagrange multiplier A(x) for the ensuing finite 
dimensional problem exists. Moreover our numerical Lagrange multipliers stayed stable 
under refining the mesh, suggesting that f-representability holds. Establishing this rigorously 
is an interesting open problem. 

Next we describe our algorithm for solving (6.7). We drop the spin variables for sim¬ 
plicity; extension to the spin-dependent case is straightforward. Equation (6.7) looks like 
an eigenvalue problem, but the “eigenvalue” J2i=i ^( x i ) depends on a function on £7, and 
moreover, the “eigenfunction” ^ has to satisfy some nonlinear marginal constraints. Due 
to these difficulties, there is no simple way for us to solve this problem directly. If we look 
at the equation (6.7) the other way around by assuming that X(x) = A(x) with some given 
function A, then the problem is reduced to the following generalized eigenvalue problem: 
Find (i e M and i7 1 (D Ar ), such that 


a . 

--A + 
2 


N 


E c(\ Xi - Xj \) *(* 1 , = y X(xi) ...,x N ) (6.8) 


1 <i<j<N 


\i=l 


with ||^||l 2 = 1 and fi being the lowest eigenvalue. We thus obtain an eigenfunction 'k with 
corresponding single-particle density 


p(x) — N 


f 


V(x,x 2 , 


,x N ) 


2 

dx 2 • • • dxjsf. 


We denote the above process (from A to p) by T, that is, p 
equivalent to the nonlinear problem 


^(A). We have that (6.7) is 


Po = J r ( A). 

We resort to the following Newton algorithm for solving this nonlinear problem. 


(6.9) 


Algorithm 6.1. Newton algorithm for solving (6.9) 

1. Fix e > 0. Let k = 1 and initialize X\ ^ 0. 


2. Solve (6.8) with A = A& to obtain (pk,^k)- 

3. Let \' k = p k \ k and p k {x) = N f QN -! \'& k (x,x 2 , ■ ■ ■ ,x N )\ 2 dx 2 • • -dx N . 

4 . If \\po — Pk || < £, get ^ and goto 5; else, let 


Xk+1 ~ X ' k + ( M 


-l 


(pO - Pk)- 


( 6 . 10 ) 


A=A< 


Take k = k + 1 and goto 2. 

5. Calculate the pair density p 2 (x,y) = (^) f QN -2 \^k( x , y, £ 3 , • • • , xj\f)\ 2 dxs • • • dx^. 

Note that the operator ^ | A _ Afc ^ in ( |6.1Q| ) can not be obtained explicitly, an approx¬ 
imation for it has to be made. We abbreviate • • • ,xn) = ^k{ x i) and obtain by 

the chain rule that 


d_T_ 

d\ 


A—A k 


dp k d^k dA k 
dV k ' dA k ' d \ k ' 


( 6 . 11 ) 


The first and third factors on the right-hand side of (6.11) can be obtained explicitly. To 
calculate the second term, we observe that 


H^k = Afe’Jfe with H = --A 


E 


- •oD- 


( 6 . 12 ) 


1 <i<j<N 
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By differentiating (6.12) with respect to and ignoring the A^-dependence of 'I ' k on the 
right-hand side, we can obtain the approximation 


dA k 


n-^k. 



obtain the true pair densities for different coupling constants. We perform all our following 
computations in double precision arithmetic on a PC with 16GB RAM using Matlab. 


6.1 Bosons 

To elucidate pure correlation effects undiluted by exchange, we first neglect the spin variables 
and assume that T is symmetric, that is, we assume that the particles under consideration 
are bosons. 

Let T be a partition of Q = [—L, L] (L = 5) with equally spaced nodes a 1 < a 2 < • • • < 
a 171 . Denote by Xj( x ) the piecewise linear function with value 1 at node a- 7 and 0 otherwise. 
Then the functions 


{Xj(x) : j = !,••• ,m} 


form a linear finite element basis set on D, which gives a discretization for the single-particle 
space. Denote the finite dimensional space spanjxj : j — 1, • • • , m} by V m . 

Since the wavefunction T in (6.4) is a function on Q N , we shall generate a basis set in 


TV-particle space by taking tensor products of the {Xj} : 


N 


*S( X ) = II*.(*‘) with J = Jn) e {!,"• ,m} N . 


(6.13) 


1 = 1 


Note that the number of degrees of freedom for this basis set is m N . We denote by B ^ the 
TV-boson space spanned by the basis functions {^j}* 


With the above discretization, we have the following variational formulation of (6.7): 
Find A G V™ and T G B 1 ! such that 


N 


-(W,Vv)+ Y, «\xi-x j \)^,v) = Y(Kx i )^,v) V«eK 

1 <i<J<N i= 1 

p(x) = N |T(x, ^ 2 , • • • , x^)\ 2 dx 2 • • • dxjsi with x = a 1 , a 2 , 

l Jqn-1 


N 


(6.14) 


The second line of ( |6.14| ) is a discretization of the marginal constraint, which is only imposed 
on the nodes of T. Within this discretization, p — pk in Algorithm |6.1| is calculated as a 
vector {p(ai) — pk(ci^)}'JL 1 on the nodes. 


For the homogeneous density ( |6.1| ) and the inhomogeneous density (6.2) with TV = 2, 3, 

with 5 = 10 -4 and 


and 4 particles, we compute their pair densities by using Algorithm 6.1 


m = 40 for TV = 2, 3, m = 32 for TV = 4. The results for different values of a are presented 
in Figure [63] and [R6] When a = 0, the electrons are strictly correlated to each other: the 
position of one electron fixes all positions of the other electrons, and the pair densities are 
given by ( |6.3| ), with support supp(pf C£; [p]) = {(^,T^(x)) : x G fi, i = 2,TV}. To visualize 
this limiting pair density, we plot, above each curve {(x,Ti(x)) : x G D}, the prefactor 
p{pc) /y 1 + (T/) 2 of the normalized line element ds along the curve. 

We observe that when a is small (e.g., a = 0.1), the pair densities are highly localized 
as 2(TV — 1) ridges around supp(pf C£; )- As a increases, the pair densities are smoothed out 
gradually. The 2 (TV — 1) ridges are still visible when a = 1 but merge with each other when 
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N = 2 


N = 3 


TV = 4 




-5 -5 -5 -5 -5 -5 


Figure 6.5: (bosons) Pair densities with p(x) = ^ (eqs. (5.1), (6.6), but with spinless 
symmetric 4/). Top row: SCE/optimal transport (based on exact results [57]). 
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N = 2 


N = 3 


TV = 4 




Figure 6.6: (bosons) Pair densities with p(x)m^j- (l+cos(^x)) (eqs. ( |5.1| ), (6.5), but with 
spinless symmetric 4/). Top row: SCE/optimal transport (based on exact results [37]). 
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a =10 and 100. The profiles of the pair densities strongly reflect the number of particles 
(particularly when a is small), a phenomenon that is missed by the standard DFT models. 
When a equals 100, the pair densities are ver y clo se to the statistically independent function 


\ (l — p(x)p(y) predicted in Proposition 5.2 In fact, the behavior of the pair densities 


as a goes from 0 to infinity can be viewed as a process in which the Coulomb holes fade 
away and the correlations are smoothed out towards statistical independence. 

Moreover, we plot the Lagrange multipliers X(x) for systems with 4 bosons and different 
values of a in Figure 6.7 We have mentioned that —X(x) can be viewed as the external 
potential that has p as the ground state density. Therefore, shifting X(x) by an additive 
constant makes no difference, and we can use an appropriate shift to allow better comparisons 
in the picture. When a = 0, the SCE Lagrange multiplier can be calculated according to 
the formulae in urn]. When a is small, the potentials are actually quite close to the SCE 
case. As a increases, the potentials converge to constant functions for homogeneous systems 
and become steeper and steeper for inhomogeneous systems to cancel the kinetic energy and 
constrain the particles. 




Figure 6.7: The Lagrange multipliers X(x) for systems with 4 bosons. Left: p(x) = 
Right: p(x) = ^ (l + cos(Jx)). 


6.2 Fermions 

Let us now come back to the fermions with spin variables and antisymmetry constraint in 
An- For simplicity, we use the notations j" and for spin up and spin down, respectively. 
With the same partition T of £4 as that in Section |6.1[ the single-particle basis set becomes 

{ Xj,s(x, a) : j = !,-■■ ,m, s =f, |}, 


where x G fl, a G and Xy s (x,cr) = Xj{ x ) if s = a and 0 otherwise. The classical 

product (6.13) needs to be replaced by the Slater determinant of the N one-body 

basis functions x il?Sl ,...,x iiV?SAr , with j = (ji, • * • J N ) e {1, • • • ,m} N and s = (si,--- ,%) G 
The number of degrees of freedom is ( 2 ^). We denote by the Wfermion space 
spanned by the basis functions {V’ys}- 

The corresponding variational formulation of ( |6.7| ) now reads as follows: Find A G 
and T e ViY such that 


N 


-(Vtf,V«)+ ]T (c(|x i -^|)«', v ) = ^(A(^)«', v ) WveV, 

1<z<j<7V 


N 


i=1 


p(x) = N [ 

jv'S!"- 1 




<7i, x 2 , cr 2 ,■■■ , x N , cr N )\ 2 dx 2 ■ ■ ■ dx N 
with x = a 1 , a 2 , • • • 


(6.15) 
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Using Algorithm |6 .1 1 with 5 = 10 4 and m = 40 for N = 2, 3, m = 32 for AT = 4, we calculate 
A and T for homogeneous and inhomogeneous electron densities given by (l6.ll) and (6.2). 


The ground state pair densities are depicted in Figures 6.8 and 6.9 As an illustration of 
the efficiency and stability of Algorithm 6.1, we present a convergence curve of ||po — Pk\\ h 1 
Figure |6.10| 

First of all, when N = 2, we have the same pair densities as those of bosons. This is 
easy to understand since for two-particle systems, the two spin variables are always paired 
up, and the antisymmetry constraint does not affect the spatial variables. 

We also find similar pair densities for bosons and fermions when a is small (e.g., a = 0.1). 
In this case, the particles are strongly correlated to each other for both bosons and fermions, 
and are always localized in different regions of space that have very little overlap. Therefore, 
the pair densities are almost independent of the choice of the spin variables: both the 
symmetric and antisymmetric choice give very similar spatial distributions, and the particles 
do not sense very much whether they are fermions or bosons. From the pictures, we can also 
draw some similar conclusions as those for bosons: When a is small, the particle number 
can be recovered by counting the number of ridges of the pair densities. As a increases, the 
2 (TV — 1) ridges merge together. 

A significant difference between bosons and fermions is that, when a goes towards infinity, 
the pair densities of fermions do not become statistically independent if AT > 2, but are 
depleted near the diagonal x = y, a phenomenon known as “exchange holes” (see e.g. ED). 
As a increases, the effects of Coulomb repulsion get weaker and weaker and the Coulomb 
holes are fading out, whilst the exchange holes take over. For comparison, the theoretical 
p2 as a —)> oc (for homogeneous p with AT = 4) is plotted in Figure 6.11 It corresponds to 


a Hund’s rule selection from the degenerate ground state of T (see Theorem 7.2), consisting 
of the orbitals 0t, 0|, It, (—l)t (hi the notation (7.16)). We observe that it is extremely 


close to the numerically computed pair density at large a (shown for a = 100 in Fig. 6.8 
bottom right panel). 

The Lagrange multipliers for systems with 4 fermions are presented in Figure |6.12| In 
comparison with those of bosons, they also converge to a constant potential as a increases for 
homogeneous systems, and have a steeper potential at the same value of a for inhomogeneous 
systems. 

From the numerical simulations in Section [6T] and [631 we conclude that the pair densities 
across the whole range of coupling constants are deformed versions of the two limit cases 
a = 0 and a = 00, with a slow and steady cross-over and without any additional effects 
appearing. The “information” in the pair densities for all a can somehow be recovered from 
just the two end values a = 0 and a = 00. By contrast, none of the end-value pair densities 
gives useful information about what happens at the other end. This lends theoretical support 
to the idea in [391 of two-end interpolation functionals. It should be very interesting to try 
to relate the specific functional proposed there to an underlying pair density model and 
compare to a theoretical adiabatic connection curve. 

Let us also emphasize the strong pair density localization without single-particle local¬ 
ization and the strong AT- dependence. The latter is missed completely by the local density 
approximation (LDA), which is based on uniform electron gas theory (AT = 00). As regards 
the former effect, it is not clear (at least to the authors) to what extent it is accounted for 
by any of the models used in practice. For homogeneous p and large AT, the true pair den¬ 
sity profile is captured implicitly through use of the LDA correlation energy; but we do not 
know what happens implicitly to the pair density when applying, say, the LDA or gradient 
corrections or a fraction of exact exchange to a typical inhomogeneous p. See Section [8] for 
further discussion. 

To end this section, we summarize some of the characteristics of the pair densities in 
Table |TJ 
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N = 2 


N = 3 


TV = 4 




-5 -5 -5 - 5 -5 - 5 



Figure 6.8: (fermions) Pair densities with p(pc) = ^ (eqs. (|5.1| ), ( |6.6[ )). Top row: 
SCE/optimal transport (based on exact results [37] . see also jT2l 111] ) . 
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N = 2 


N = 3 


TV = 4 




Figure 6.9: (fermions) Pair densities with p(pc) = ^ (l + cos(^x)) (eqs. ( |5.1| ), ( |6.5| )). Top 
row: SCE/optimal transport (based on exact results [37] . see also mm)- 
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Figure 6.10: Convergence curve of \\pk — po\\ 
of Algorithm |6.1 1 for N = 2 and a = 0.1. 


Figure 6.11: Theoretical pair density as a 
oo, coming from the spin-polarized Slater de¬ 
terminant |0t OJ, It (—l)t) (Theorem 


7.2). 




Figure 6.12: The Lagrange multipliers \{x) for systems with 4 fermions. Left: p(x) = 
Right: p{pc) = ^ (l + cos(Jx)). 


7 Rigorous asymptotic results 


The following asymptotic results in ID support our numerical findings, and were used to 
test the correctness of our code. Results of this type are well-known in the physics literature 
(except perhaps those on “selection rules” which emerge in the non-interacting limit in case 
of orbital degeneracies) and the novelty consists only in providing rigorous proofs. The 
reader is reminded that on the rigorous level very little is known about exact DFT and even 
basic issues as raised in m such as continuity of the HK functional remain open. 

Recall from Section [ 5 ] the scaled density-to-pair-density map p2, a [p\ = D a -ip 2 [D a p\, 
where p 2 is the original density-to-pair-density map. 


Theorem 7.1. (Small a limit, ID systems) Let p be any single-particle density on M be¬ 
longing to the class IZn (see ( |2.9 | ) ), N > 2. Assume that p > 0 in some finite or infinite 
interval (a,b), and p — 0 outside. Let Ti,..,T/v be the following optimal transport maps 
found in f37^ and justified rigorously in m for N = 2 and in m for general N: let 
do = a < di < ... < djsi -1 < d n = b be the partition of (a,b ) into N sub-intervals of equal 


mass, i.e. 



1 (i = l,...,A0, 


and let X 2 be the unique p-preserving map which monotonically maps each interval [di-i,di\ 
(i=l,...,N-l) to the next interval 1 ] and the last interval [d/v-i?d/v] to the first, [do,di\. 

Let T\(x) = x, and let Tj, j = 3,..,N, be the (j-l)-fold composition of with itself. (See 
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pair densities 

bosons 

fermions 

o = 0 

SCE 

SCE 

0 = 00 

statistical independence 

single slater determinant 

Coulomb holes 

yes 

fade out as o increases 

yes 

fade out as o increases 

exchange holes 

no 

yes 

fade out as o decreases 

7V-dependence 

yes 

yes 


Table 1: Summary of the characteristics of pair densities. 


Figures 6.3\ 6.J\ ) Then 


1 N 

^P2Ap\ = ~ 2 Y. 

3 =2 


p(x) 


1 + T'Ax) 2 


-ds 


y=Tj (x) 


the limit being in the sense of weak* convergence of Radon measures. 


Proof Let be a minimizer of the variational problem in (5.1), and let pN,a( x u --, x n) = 
^2 Sl s N ez 2 \^ot( x i, s u x Ni $n) | 2 - Since the pN,a have marginal p, they are a tight family 
of probability measures (to show this one proceeds analogously to the proof of a similar result 
in the appendix of [29]) and hence possess a subsequence (see [O]), again denoted pN,a, 
converging weak* to a probability measure p^v,* as a -A 0. By standard arguments p ^,* 
has one-body marginal p. Moreover, by dropping the kinetic energy from (5.1) and using 
the lower semicontinuity of the interaction energy under weak* convergence, and letting 
Vee = T,i<n c Ai ~ x j), 


lim inf (aThEl + V ee m) = lim (aT[^ 0 
a -^0 L J L a —>-0 V L 


+ 


/ V ee pN,a) > / V ee dpN,*- 

JRN / J r n 


On the other hand, as proved in [12] the left hand side equals min PAr f V ee dpN, the minimum 
being over symmetric probability measures on with marginal p. It follows that p^/-,* 
is a minimizer of the latter problem. By the results of m as made rigorous in m, the 
minimizer of the latter problem is unique and given by ( 4.15| ), with the above explicit maps 
Ti,..,T/v- The uniqueness implies that the whole sequence pN,a converges weak* to Pat,*- 
Next, this latter convergence implies weak* convergence of the associated two-body density 
pJf Q to the pair density (^) f pN,*dx 2 .--dxN of pw,*- The assertion now follows from our 
result (|6.3|) . □ 


Here and below, we denote the eigenfunctions of the Laplacian on [—L, L\ with periodic 
boundary conditions by 


I k){x) := 


Ak^-x 


^/2 L 


(he Z) 


(7.16) 


and the associated spin-orbitals \k)(x)S^(s) and k)(x)Si(s) by \k t), \k i). 


Theorem 7.2. (Large a limit, homogeneous ID systems) Let p(x) = p = N/(2L) be the 
homogeneous density on [— L,L\, and let P 2 ,a[p\ = D a -ip 2 {D a p\ be the scaled density-to-pair- 
density map for periodic boundary conditions on [—L,L\ ( ( |5.1| ) with An given by ( |6.6 | ) ). Let 
T be the Slater determinant built from the first N orbitals 0i,..., <pN of the (partially spin- 
polarized) sequence |0 t), |0 |1 t), |(-1) t), I 1 I), |(~1) I); I 2 t>, l(~ 2 ) t), I 2 l), K” 2 ) I); 
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.... Then, letting z = j^(x — y), 


lim p 2 ,a[p](x,y) = P2(x,y) 

ot-t-oo 


1 _ 2 1 sin 2 (^z) 

2 P ~ (2i) 2 sin 2 (^) ’ 

1 _2 11 sin 2 (^2:)+sin 2 (^±±2:) 

' 2 P “ 2 (2L) 2 sin 2 (iz) 

1 _ 2 1 1 sin 2 (^,z) + sin 2 (^±2,z) 

X “ 2 (2L) 2 sin 2 (± 2 ) 


X = 2 mod 4 
X = 1 or 3 mod 4 
TV = 0 mod 4, 


t/m limit being in the sense of strong convergence in L 1 ([—L, L] 2 ). 

Proof. We first ignore the constraint 4T i—)► p. Let Xo be the ground state of T = — on 
A at, let Po be the orthogonal projector from L 2 onto Xo, let X^ be the lowest eigenspace of 
PoV ee Po within Xq (note that Xq = Xq if Xq is one-dimensional), and let S ' 0 = {T E Xq : 
4/ i—^ p}. By degenerate first-order perturbation theory, together with the fact that by the 
explicit description below S ' 0 is nonempty, 


lim {4/ G An | 4/ minimizes T + 

a—>- oo 


i s/to p} C S' 0 , 


(7.17) 


the limit being in the sense of strong L 2 convergence. It follows that the set of pair densities 
p 2 ,a[p\ satisfies lim^^oo p 2 ,a[p\ ^ {p 2 : 4/ £ Pq}, th e limit being in the sense of strong 

L 1 convergence (note that the map 4/ 1 —p 2 is continuous from P 2 (([—L,L] x Z 2 ) N ) to 
P 1 ([—L, L] 2 )). To complete the proof of the theorem, we need to understand S ' 0 explicitly. 
The ground state Xo of — on An is given by 

Span |0 2 , l 2 , (-1) 2 ,., X 2 , (-X) 2 ) if X = 2 mod4, K = (7.18) 


and by 


Span{|0 2 , l 2 , (—l) 2 ,(X-l) 2 , (-(X-l)) 2 ,ai, ..,a d ) : a 1 ,..,a d = any d 

orbitals from —K j", —K otherwise, (7.19) 

where the notation k 2 means that the orbitals \kf) and \k\) are both present in the Slater 
determinant and d and K are as follows: d = 3 and K = (X — l)/4 if X = 1 mod 4; d = 2 
and K = X/4 if X = 0 mod 4; and d = 1 and K = (X + l)/4 if X = 3 mod 4. For X^O 
mod 4, Xq = X 0 . But for X = 0 mod 4, aligning the two spins is favourable because it 
generates an additional exchange term. This is a manifestation of the empirical Hund’s rule. 
Thus Xq is given by the subspace of Xq with total spin S 2 = s(s + l)| s =i, 


X' = Span{|0 2 , l 2 , (-1) 2 , (K- 1 ) 2 , -(K- 1 ) 2 ,K f, -K f), 

|0 2 , l 2 , (-1) 2 , (K- 1) 2 , -(K-1) 2 ,K i, -K |), 

^(|0 2 , -{K- 1) 2 , K t, (-K) i) + |0 2 , -(K- 1) 2 , K 4., ( -K ) t))} 
if X = 0 mod A. (7.20) 

The three states above are the canonical basis states with S 3 = 1, —1, and 0. 

We now take into account the constraint T p, and determine S' 0 . For even X, S ' 0 is 
the sphere of unit vectors in Xq. For odd X, we claim that 

S'o = {a*i + P*2 + 7*3 + <5*4 : H 2 + |/?| 2 + 111 2 + H 2 = 1, (jf) ■ (j) = 0}, (7.21) 

where for X = 3 mod 4 the 4/i,.., 4 /4 correspond to the four choice of a\ in ( |7.19| ) in the 
listed order, and for X = 1 mod 4 they correspond to the four choices K ^ (—X) | X 
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K t (-K) t K i, K i ( -K ) J. (-K) f, and K t (-if) f (-if) J, of m, a 2 , a 3 . For, say, the 
latter N’s, the constraint in (7.21) follows from the fact that 


p a* 1 +."+i* 4 ( x ) = consi+ | Q;e iKf a:+7e -iiff a: |2 + | /?e iKf a: + (5e -iKf a: |2 

= const + 2 Re (cry + /35) cos( 2 ATy x) — 2 Im(cry + /35) sin( 2 IT yx) 

±j 1 j 

and the linear independence of the three functions cos, sin, and 1. Finally, for each of the 
four cases of APs, a tedious calculation gives the corresponding pair densities, as well as the 
fact that these are independent of the coefficients of the wavefunctions in S' 0 . □ 

We find the uniqueness of the limiting /^’s despite degeneracy of the limiting ground 
state wavefunctions remarkable. 


8 An ansatz for homogeneous systems 


Based on the above numerical and asymptotic results, we shall now design a simple ansatz 
for the pair density of homogeneous systems which is accurate across the whole range of 
coupling constants a. 

If we look at the pair density graphs for homogeneous systems from a specific angle (see 
Figure 8.1 for example), we can observe that they are almost uniform functions of x — y. 
This together with the peaks on the graphs of the transport maps Ti suggests an ansatz of 




Figure 8.1: Rotating the pair density of a homogeneous system with 4 fermions with a = 1. 
Left: view from angle (-35,50). Right: view from angle (-45,0). 


the form 


P2(x,y) 



T(di{x,y)) 


where di(x,y ) = min|(x,y) - (x , ,T i (x , ))|- 

x' 


( 8 . 1 ) 


Here c n is a normalization constant and T is some shape function. Note that, due to the 
explicit form of the the above p 2 depends only on x — y. A general formal asymptotic 
expansion at small a in the physics literature [ 20 ] or alternatively, in our special case, an 
elementary calculation detailed below suggests to take T to be a Gaussian. Thus we make 
the ansatz 


P2(x,y) « G^ os (x,y) = c n ^XP x p(~ ^^ ( 8 ‘ 2 ) 

where the parameter is allowed to depend on the coupling constant a and the particle 
number N. To obtain <7 we minimize the id-error \\p 2 — G^ os ||^i(^ 2 y where p 2 is the correct 
pair density as computed in Section [ 6 ] 

See Table [ 2 ] for the optimal parameters c as well as the error (in different norms, calculated 
by using the finite element discretizations used in Section^ between the correct pair densities 


25 














and the ansatz (8.2). We present some cross sections (on x = —y) of the pair densities 
and our ansatz in Figure 8.2 It appears that the ansatz ( |8.2| ) provides quite an accurate 


approximation. Note that the ansatz (8.2) is accurate at the two limits (by taking <; = 0 at 
a = 0) and (s = oo at a = oo), and we can observe from Table [2] that the approximations 
are better in the regimes where a is very small or large. 

Finally, we give the promised elementary argument which lends theoretical support to 
our Gaussian ansatz. For a = 0, c(r ) = 1/r, and, say, N = 2, the Lagrange multiplier in eq. 
(6.7) is known exactly and equals \(x) = \x\/L 2 . Hence the total potential in (6.7) is 


V(x,y) = 


1 


\x-y\ 


FI 

L 2 


L 2 ' 


This potential is minimal on graph = {x — y = d= L}. For nonzero but small <a, the 
ground state should still be localized near graph T2, and hence we may replace V(x,y) by 
its second order Taylor polynomial at the nearest point to (x,y) on graph T2. This Taylor 
approximation is easily calculated to be 

f>/~„\_ 2 , d 2 (x,y) 2 _ 2 | min {(x - y - L) 2 , (x - y + L) 2 } 

V(x,y) l+ l 3 L + L 3 

Eq. (6.14) with this potential is solved exactly by a Gaussian of form e~ d ^ x,v ' > / const j except 
on the diagonal x = y, where the Gaussian and the exact solution should both be small 
and hence close to each other. This suggests that eq. ( |8.2| ) (with N = 2) is a good global 
approximation to the pair density. Giving a rigorous version of this argument is an interesting 
open problem. 


N 

a 

optimal <^ 2 

\\P2~G*»\\ L i 

\\P2 - G*° s \\ L 2 

v ee [p 2 ] 

V ee [p 2 ] - V ee [G b ™} 

2 

0.1 

1.21 

0.0563 

0.01453 

0.218 

-0.00742 

0.3 

1.79 

0.0773 

0.01126 

0.243 

-0.00876 

1 

2.70 

0.0604 

0.00782 

0.277 

-0.01147 

3 

3.78 

0.0472 

0.00662 

0.339 

-0.00662 

10 

6.86 

0.0026 

0.00358 

0.420 

-0.00605 

100 

52.1 

0.0013 

0.00097 

0.667 

-0.00107 

3 

0.1 

0.67 

0.0922 

0.01850 

0.814 

-0.01841 

0.3 

1.01 

0.1146 

0.02284 

0.861 

-0.02098 

1 

1.42 

0.1571 

0.03816 

0.932 

-0.03367 

3 

2.19 

0.1867 

0.02122 

1.170 

-0.01880 

10 

7.02 

0.1292 

0.01245 

1.438 

-0.01138 

100 

64.0 

0.0388 

0.00469 

2.009 

-0.00738 

4 

0.1 

0.90 

0.2408 

0.03570 

1.899 

-0.02522 

0.3 

1.42 

0.2885 

0.03687 

2.092 

-0.04871 

1 

1.93 

0.2809 

0.05231 

2.162 

-0.05275 

3 

9.48 

0.3788 

0.04266 

2.773 

-0.05477 

10 

32.1 

0.1024 

0.02387 

3.234 

-0.02013 

100 

232.0 

0.0542 

0.00681 

3.751 

-0.00045 


Table 2: Approximations of the pair densities of homogeneous systems (for bosons). The 
Coulomb energy in the last two columns is defined by V ee [p 2 \ = f f P 2 (x,y)c(jx — y\)dxdy. 

For fermions, to capture the asymptotic emergence of exact exchange as a -> oo we make 
the ansatz 

G l ^r,(x,y) = c n G* os (x,y)^ (p(x)p(y) - r]r(x,y )), (8.3) 

where c n is a normalization constant, 77 G [0,1] is a parameter (allowed to depend on N 
and a), and r is the exchange term from (|4.9|). The freedom of varying p allows a seamless 
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Figure 8.2: The cross sections (on x — —y ) of the pair densities and their approximations G^ os for 
homogeneous electrons. 


crossover between the SCE pair density (rj = 0, c = 0) and the exact-exchange pair density 
(rj = 1, s = oo). The ansatz (8.3) is not the only way to achieve this, but it is perhaps 
the simplest. Note that, unlike in B3LYP |3], exchange is mixed in multiplicatively, not 
additively. Numerically, we obtain 77 by minimizing the L 1 -error \\p 2 — G^Wl 1 ^) ( while 
keeping, for simplicity, the bosonic values of ^). The results in Table [ 3 ] and Figure [873] show 
that (8.3) is a good approximation for fermions. In particular, Figure |843| (which concerns 
the case N = 4 and different values of a) shows that the transition from 6 (= 2 (TV — 1 )) 
SCE ridges to 4 exact-exchange ridges is correctly captured. The ansatz ( |8.3| ) is accurate at 
the two limits a = 0 and a = 00 , and the approximations are indeed better in the regimes 
where a is very small or large, as we can see from Table [ 3 ] Moreover, we observe that the 
errors for fermions are larger than those for bosons, which may be caused by the complicated 
interplay of Coulomb and exchange holes. 


N 

a 

optimal r] 


\\P2-G^\\ L 2 

VM 

VM - Vee[G*%\ 

3 

0.1 

0 

0.1375 

0.02267 

0.814 

-0.01694 

0.3 

0 

0.2030 

0.03601 

1.416 

-0.01173 

1 

0.01 

0.1919 

0.02706 

0.926 

-0.04983 

3 

0.02 

0.2217 

0.02681 

1.094 

-0.04810 

10 

0.27 

0.1792 

0.02139 

1.345 

-0.04138 

100 

0.92 

0.0264 

0.00324 

1.676 

-0.00761 

4 

0.1 

0 

0.3262 

0.03189 

1.898 

-0.04645 

0.3 

0 

0.3222 

0.03991 

2.061 

-0.07957 

1 

0.01 

0.3457 

0.04296 

2.133 

-0.07976 

3 

0.03 

0.3539 

0.03961 

2.675 

-0.08353 

10 

0.47 

0.1449 

0.01752 

2.971 

-0.06630 

100 

0.95 

0.0223 

0.00306 

3.134 

-0.00155 


Table 3: Approximations of the pair densities of homogeneous systems (for fermions). 


9 Conclusions 

In this paper we studied the exact density-to-pair-density map in density functional theory. 
In the absence of any previous numerical simulations of this map, we computed it here 
for typical one-dimensional families of densities obtained by scaling. This is the same as 
computing the map along the (two-sided) adiabatic connection from the non-interacting 
limit to the strictly correlated limit. We observed a slow and nontrivial cross-over between 
the endpoint profiles, which are given by exact exchange respectively by SCE correlations (or 
mathematically: by first-order perturbation theory respectively by optimal transport with 
Coulomb cost). The cross-over, while smooth, is very far from a linear interpolation and 
involves multiple lengthscales. 
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Figure 8.3: The cross sections (on x — —y ) of the pair densities and their approximations for 
homogeneous systems with 4 electrons. 


This study gives us a deeper insight into the details of electron correlations, and may 
further lead to novel models for the pair density (and hence the interaction energy). As 
a fist step, we constructed an ansatz for pair densities of homogeneous systems in one di¬ 
mension which is exact in the weak and the strong interaction limit and has been shown to 
remain accurate in the whole intermediate regime. The ansatz itself is readily generalized 
to inhomogeneous three-dimensional systems, but for such systems we have not yet tested 
its accuracy in the intermediate regime, nor do we know how to pick the correct parameter 
values just from the one-body density. We hope to come back to these issues in future work. 
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